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Abstract 



In this note we propose and analyze novel implicit-explicit methods based on second 
order strong stability preserving multistep time discretizations. Several schemes are 
developed, and a linear stability analysis is performed to study their properties 
with respect to the implicit and explicit eigenvalues. One of the proposed schemes 
is found to have very good stability properties, with implicit A-stability for the 
entire explicit stability domain. The properties of the other proposed schemes are 
comparable to those of traditional methods found in the literature. 
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1 Introduction 

Implicit-explicit schemes (IMEX) are methods for the solution of time-dependent 
differential equations in the form 



where / and g are terms with different character such that one, say g, man- 
dates implicit treatment whereas y = f can be solved efficiently by an explicit 
method. Such systems often arise from spatial (semi-) discretization of time 
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dependent PDEs by the method-of-lines - a prime example being discretiza- 
tions of transport equations that may contain different terms accounting for 
advection, diffusion, and reaction. 

It is beyond the scope of this note to give a complete description of IMEX 
methods; for a systematic discussion of multistep IMEX schemes we refer to 
Ascher et al. [2]. Frank et al. [3] analyzed the stability of multistep IMEX 
schemes, and showed that stable implicit and explicit integrators do not nec- 
essarily lead to a stable imphcit-explicit method. Ascher et al. [1] have also 
developed IMEX schemes based on multistage Runge-Kutta integrators. 

The Strong Stability Preserving (SSP) property is a generalization of the To- 
tal Variation Diminishing (TVD) property for hyperbolic problems, and it is 
a guarantee of non-linear stability of the scheme. Traditional explicit multi- 
step schemes, such as the Adams-Bashforth methods, do not necessarily have 
the SSP property [4]. Note however that recent work - by Hundsdorfer and 
Jaffre [6] and by Hundsdorfer, Ruuth, and Spiteri [8] - have shown that the 
traditional multistep methods can be made monotonicity-preserving, e.g. with 
respect to the Total Variation semi-norm, if a suitable starting procedure is 
employed. 

In this note we discuss the linear stability of novel IMEX schemes based on 
second order accurate SSP multistep time discretization. Although SSP is a 
non-linear stability property, it is nevertheless appropriate to study the linear 
stability of the proposed schemes, since IMEX combinations do not necessarily 
inherit the stability properties of the component schemes [3]. A full analysis 
of the nonlinear stabilities of the proposed schemes is beyond the scope of this 
note. 

We will show that second order IMEX schemes based on the explicit multi- 
step SSP discretizations proposed by Shu [11] have stability properties that 
are comparable to the Crank-Nicholson/Adams-Bashforth and BDF2 IMEX 
methods. 



2 Stability analysis 

A general implicit-explicit k-step method for (1) can be written in the form 

k k k 

X] (^iVn+l-i = XI ^ifn+l-i + X] (^i9n+l-i, (2) 
i=0 i=l i=0 

where /„ = t„) and Qn = g{yn: tn)- 
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The stability analysis below is based on the scalar test equation 



y^Xy + ny, (3) 

where A and fi arc complex constants that represent the eigenvalues of the 
explicit and implicit operators, respectively. 

By applying the multistep method (2) to the scalar test equation (3) and 
looking for solutions in the form y = C", we obtain the characteristic equation 



1=0 1=1 i=0 

For convenience, we follow Prank et al. [3] in working on the transformed form 
of the equation 

A(z) - \B{z) - ixC{z) = 0, (5) 

where z = 1/C, and the polynomials A, B, and C correspond to the time 
derivative, explicit, and implicit operators, respectively, and are given by 

A{z) = j2ai^\ B{z) = j2biz\ C{z) = j2ciz\ (6) 

1=0 1=0 i=0 

Note that, since we are working in the transformed variable z, the scheme (2) 
is stable if all the roots of the characteristic equation (5) are in the exterior 
of the unit disk, |2;| > 1 (with strict inequality if 2; is a multiple root). 

To investigate the stabihty properties of the combined IMEX scheme we will 
study the image of the exterior of the unit circle under the mapping 

, , A{z) - \B{z) 
C[z) ■ 

The boundary of the stability domain, S, for the explicit method can be readily 
determined from 

dS^{X: X{0) = A{e^^)/B{e^^), -tt < 9 < tt} . 

Similarly, we can find the stability region, V, of the implicit method from 

&D^{fx: n{e) = A{e^^)/C{e^^), -tt < ^ < tt} . 
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3 IMEX methods based on second order multistep SSP schemes 



3.1 SSP time discretizations 



The first order accurate forward Euler explicit time integration is strongly 
stable for some norm, || • ||, if 

hn + ^tfnW < IbnII, 

under a suitable time step At < Ato, where Ato is the largest allowable time 
step for stability. Higher order time discretizations are strong stability pre- 
serving (SSP) if they retain this property, possibly for a shorter time step 
At < cAto- The constant c < 1 is, by convention, called the CFL coeffi- 
cient [4], and should not be confused with the Courant number of hyperbolic 
discretizations. The SSP time discretizations were originally developed by Shu 
[11] for the Total Variation norm to maintain positivity and monotonicity in 
the solution of hyperbolic equations. 

Shu [11] proposed a family of exphcit SSP multistep schemes. We will restrict 
the discussion to second order schemes with positive coefficients. We thus 
consider the three- and four-step second order schemes: 

. [3] _ - 3yn - yn-2 _ . . 

. [4] _ Qyn+1 - Syn - yn-3 _ r x 

yM /n, (8b) 

where denotes the fc-step discretized time derivative for the SSP scheme. 
These methods are strong stability preserving with CFL coefficients c— 1/2 
and c = 2/3forA; = 3 and k — A, respectively. The polynomials (6) that 
characterize these schemes are 

A^^\z) = (4- 3z — z )/6, (9a) 

^W(^) = (9_8^-^4)/i2, (9b) 

and 

B{z) = z. (10) 

We show the explicit stability domain, S, for the two schemes (8) in Fig. 1. 
Notice that both schemes have an appreciable part of the stability region 
close to the imaginary axis, something that is advantageous for the solution 
of hyperbolic terms. 
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Fig. 1. Explicit stability domain, S, for the two SSP time discretizations (8). 
3.2 Construction of implicit integrators 

To develop IMEX methods based on the multistep SSP schemes, we first 
construct implicit schemes by replacing the term /(?/") in (8) above by linear 
combinations of the implicit operator evaluated at the available time levels 



such that the resulting formula for calculating y""^^ is second order accurate at 
time level n. We must obviously require cq ^ for the scheme to be implicit, 
and furthermore that J2cj — 1- 

Evaluating the order conditions we obtain an under-determined system that 
allows infinitely many solutions with several free parameters. We will not go 
into all the possibilities here. Instead we focus on the straightforward approx- 
imations 



k 



j=0 




(11a) 




(1 - P)gn+i + 2Pg^ + (1 - /3)^„_i 



(lib) 



where /3 < 1 is a non-negative algorithmic parameter. 



The polynomials (6) that correspond to the right-hand-sides in (11) are then 

Ci(z) = (2 + z=^)/3, (12a) 
C2{z) = ((1 -(5) + 2(5z + (1 - (5)z^) /2. (12b) 



An implicit method has A-stability if it is stable for all eigenvalues, /i, in 
the right half plane, Re(/x) < 0, and it has 74(Q;)-stability if it is stable for 
eigenvalues in a wedge-shaped domain about the negative real axis with half- 
angle a. 

In Fig. 2 we show the stability regions of the A: = 3 versions of the implicit 
integrators (11). Note that the scheme (11a) appears to be A-stable, whereas 
the other variant (lib) is A(Q;)-stable. We can show that these stability prop- 
erties carry over to the k — A versions of the schemes. That is, (11a) is A-stable 
and (lib) is A(Q;)-stable, for both A; = 3 and k — A. 

Lemma 1 The implicit method (11a) is A-stahle for both k = 3 and k = 4. 



PROOF. Obviously, for a purely implicit method A = 0, and we consider the 
mapping 

A-stability requires that the image of the unit circle under this mapping is 
entirely in the right half plane, i.e. Re (^(po{e^^)^ > 0. 



(i) For A; = 3 we find 



Re {Me'')) = 



14- 12cos^-6cos2^ + 4cos3^ 



2 + e 



130 



5 - 6 COS 6' - 3 cos^ 6 + A cos^ 9 



2 + e 



136 



(13) 
(14) 



The denominator in this expression is obviously positive, we must therefore 
show that the numerator is non-negative. Let x — cos 9 and write 

f{x) = 4x3 - 3^2 - 6x + 5 = (4x - 5){x - if. 

It is straightforward to verify that f{x) is non- negative for —1 < a; < 1, and 
we have shown that (11a) is ^-stable for A; = 3. 

(ii) Similarly, for A; = 4, we have 



AW(ei^)M(ei^), 
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leading to 

/ / 6- llcos^ + 9cos3^-4cos^^ 
Re [(poie'")) = 



4 



2 + e- 
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Again we let x — cos 9 and write the numerator as 

f{x) = -Ax^ -9x^ -x + 6 = -(4x2 -x-6){x- if, 
which we readily verify is non- negative for — 1 < x < 1. □ 

Lemma 2 The implicit method (lib) is A{a) -stable, with a — given by 



tancco = , for /c = 3, (15a) 

(7 - 1)^ 



tanao- ^ o o ^ 3 ^ /or ^ = 4, (15b) 
2 — 37 + 7"^ 



where ^ — (3/{(3 — 1). 



PROOF. Again, we let A = 0, and we consider the mapping 

shown in Fig. 2, to determine the half-angle, ao- Note that v9o(exp(i6')) is 
singular in the limit cos^ = (3 / {(3 — 1). We only consider 13 < 1/2; since 
for P > 1/2 the scheme behaves more like an explicit method. The angle ao 
is given by the slope of the asymptote, which can be written 

tan an = hm ■. — , 

Re((^o(ei^)) 

We then have 

-4sinr-sin3r , , 

tancto = tt, for A; 3, 

4 — 3 cos 9* — cos 36'* 

-8 sin ^* - sin 4^* ^ , 

tan ctn = -7—, for k — A, 

9-8cos^* -cos4^*' 

from which the results follow. □ 



The natural choice [3 = 0, i.e. 7 = 0, gives tanag = 2 and tancto = 1 for the 
three-step and the four-step schemes respectively. In both cases, this is the 
optimal parameter value that maximizes ao- 
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Fig. 2. Stability regions, V, for the implicit integrators (11), where k = 3 and /3 = 0. 
The methods are stable outside the shaded regions. 



3.3 Stability of the IMEX schemes 



We can the construct implicit explicit methods by combining the above SSP 
exphcit schemes (8) and the implicit integrators (11). Thus we consider 



yZ = fivl + M (1 - + + (1 - /?)5(2/"^') , (16b) 



Prank et al. [3] presented two general stabihty results for multistep IMEX 
schemes. The first gives restrictions on the explicit eigenvalues in order to 
retain the full stability domain of the implicit operator, whereas the second 
gives restrictions on the implicit eigenvalues in order to retain the full stability 
domain of the explicit operator. Unfortunately, neither of these results appear 
to be applicable to the proposed schemes, and we will therefore study the 
stability properties of (16) by working directly on the mapping (px{exp{i0)) 
from (7). Remember that the IMEX schemes are stable with respect to implicit 
eigenvalues, /i, that are in the exterior of the mapping of the unit disk under 
(fx- For a method to be A-stable, the entire image of the unit circle must be 
in the right half-plane. 




(16a) 
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3.3.1 Stability of scheme (16a) 

In this section we will discuss the stability of the scheme (16a), characterized 
by the polynomials (9), (10), and (12a). For the k — ?> variant of the scheme 
we propose the following stability result: 

Conjecture 3 The scheme (16a) with k = 3 is A-stable with respect to ji for 
XeS. 

Recall that A-stability requires that the image of the unit circle under the 
mapping (px{z) must be entirely in the right half-plane for A e 5, i.e. 

$(A, z) = Re ((^A(e^^)) > 0, -n <e <7r, XeS. 

As a consequence of the maximum principle, $ will attain its minimal value 
w.r.t. A on the boundary of the explicit stability domain. We need therefore 
only consider X e dS. 

To support the assertion that (16a) with k = 3 is A-stable, we first show, 
in Fig. 3, the image of the unit circle under (px for A e dS. Furthermore, 
let A* = 74[^l(exp(i^*))/5(exp(i^*)) be an eigenvalue on the boundary of the 
explicit stability domain. Then Lpx*{cxp{i9*)) = 0, and we can investigate 
the behaviour of Lpx* close to this zero by performing a Taylor expansion. To 
leading order we obtain 

5-1-4 cos 60* 

We have neither analytically nor numerically managed to find other zeros or 
negative values for $(A, z) and the scheme does therefore indeed appear to be 
74-stable. 



For the k = 4 variant we propose: 



Conjecture 4 The scheme (16a) with k = 4 is A{a)-stable with respect to ji 
for X & S. The half-angle, a, of the wedge of stability is a ~ 0.237r. 

In Fig. 4 we show the image of the unit circle under (px for X E dS, and we 
observe that the method version appears to be A{a) stable with a 7r/4. 
Carrying out a Taylor expansion of (fx* as above, we find to leading order 

, i.,, 3 sinr -SsinSr + 2sin4r , 

^^^^ 4 sin2 3^*+ (cos3e* + 2)2 ^ ^' ^ ^ 

X / / iexx 3 6 + cosr + 3cos3r + 2cos4r , 

Im (^A ~ ^ 7 ^ ^ {d - 0*). 18b 

^^^^ " 4 sin2 3^* + (cos3^* + 2)2 ^ ' ^ ' 
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Re(z) 



Fig. 3. Image of the unit circle under the mapping ipx{z), with A G dS, for the 
schemes given by (16a) with A; = 3. The impUcit stabiUty domain of the method is 
outside the shaded region. 

We can estimate the angle for y4(Q;)-stability from the slope of the first order 
approximation: 



Numerically we find that tana ^ 0.89, and this gives a ^ 0.237r which corre- 
sponds well with the stability domain shown in Fig. 4. 



3.3.2 Stability of scheme (16b) 

We then consider the scheme (16b), characterized by the polynomials (9), (10), 
and (12b). For this scheme, we can show the following stability result: 

Lemma 5 The implicit- explicit method (16b) is A{a)-stable with respect to 
the implicit eigenvalues, /i, for P < ^ if 



tan a ~ inf 



6 + cos e* + 3 cos 3^?* + 2 cos 46* 



(19) 



sin e* -3 sin 36** + 2 sin 46** 



A 



& S'^ — S n {z : \lm{z)\ < u}, and \a\ < a, 



V 
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Re(z) 

Fig. 4. Image of the unit circle under the mapping ipx{z), with A G dS, for the 
schemes given by (16a) with k = 4. The impUcit stabihty domain of the method is 
outside the shaded region. 

where v and are given by 



for k — 3, and 



for k 



i'<(2 + 7)\/i-t73, (20a) 

(7 - 1)^ 



(2 + 7')Vl-7V3, (21a) 

tana^ — - . 21b 

2-37 + 73 ^ ^ 



PROOF. Recall from Sec. 3.2 that, in this case, the mapping ipx{exp{i9)) is 
singular in the hmit cos^ — > 7 = P/{P — 1). Taking into account a non-zero 
X — Xr + iXi & S and proceeding as in the analysis of Sec. 3.2, we have 



ior k — 3 



for k = 4. 



tancK;^ 



tan ax 




-7)yr^-3A; 
(7 - iy + 3_Xr ^ 

2 + 7^)v/r^-3A, 
73 - 37 + 2 + 3A, 
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We observe that for large |Aj|, say |Aj| > u; tanaA may become negative. 
This means that the wedge of stabihty will be rotated entirely into the upper 
or lower left quadrants. In order to have a useful stability domain that - at 
least - contains the negative real axis, tan must obviously be non-negative, 
and we must introduce an upper limit on the imaginary part of the explicit 
eigenvalues, 

<i^<{2 + 7)^1-72/3, for A; = 3 
<v< (2 + 72)^1-72/3, for A; = 4. 

We obtain a lower bound on the angle a if we choose A = ±ii/, and we have: 

(2 + 7)v/r^-3^ 

= ' 

(2 + 72)v/r^-3z/ 

tan = — - — ; — . 

2 — 37 + 7'^ 

for A; = 3 and A; = 4, respectively. □ 



Note that the natural choice = maximizes both the angle a and the upper 
bound on z/, and thus gives the largest stability domains for both the implicit 
and explicit operators. Furthermore, although the above results require that 
we use a reduced explicit stability domain, the restrictions are not prohibitive. 
The choice of the bound v can be guided by the discretization of the exphcit 
operator. As an example, we compare - in Fig. 5 - the Fourier symbol of the 
third order finite difference k = 1/3 advection scheme (see for example [7]) 
with the reduced stability domain, S^'^ for A; = 3. The Courant number in 
this particular example is cr = 0.35 which is well above the value for which the 
scheme is expected to be SSP (ctssp = 0.25). In the table below, we will use 
the value z/ = 1/3 which corresponds fairly well to cr = 0.25 for this advection 
discretization. 



4 Discussion and concluding remarks 



In this note we have proposed and analyzed the linear stability properties of 
implicit-explicit methods based on the Strong Stability Preserving multistep 
methods introduced by Shu [11]. Of the four variant methods that we consider, 
one scheme stands out from the rest. The three-step scheme (16a) with A; = 3 
retains full implicit A-stability for the entire explicit stability domain. In fact, 
we are not aware of any other second order IMEX scheme with this property. 
For the other methods, we have shown y4(a)-stability for the entire, or slightly 
reduced explicit stability domains. 
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Fig. 5. Explicit stability region S for the time discretization (8a), with the Fourier 
symbol of the third order k = 1/3 finite difference advection scheme (cr = 0.35), 
and the bound (7 = 0.5) that defines the restricted stability region S'^ in Lemma 5. 

It is instructive to compare with the traditional multistep IMEX methods 
such as the modified Crank-Nicholson/Adams-Bashforth (mCNAB) scheme 
and the IMEX BDF2 scheme. The mCNAB scheme introduces an algorithmic 
parameter, c, to the implicit integrator and can be written 

Vn+l - Z/n = l^t(3fn - fn-1 + (1 + c)gn+l + (1 - 2c)g'„ + CQn-l) ■ 

Note that c — recovers the original Crank-Nicholson scheme for the implicit 
integrator, whereas other choices of c may improve the stability properties of 
the method, see [2,3] for a discussion. The IMEX BDF2 scheme is given by 

3yn+l - ^Vn + Vn-l = 2At(2/„ - fn-1 + 5'n+l)- 

Both these schemes arc restricted to implicit 74(Q;)-stability if we require full 
explicit stability [3]. Wc summarize the implicit stability properties for all 
these methods in Table 1. Furthermore, all the methods appear to have com- 
parable non-linear stability properties. Whereas the schemes proposed in this 
work are based on strong stability preserving methods, both the Adams- 
Bashforth method and the extrapolated BDF2 method was recently shown 
to have good monotonicity-preserving properties provided suitable starting 
procedures were used [6,8]. Finally, it is difficult to assess the computational 
cost of each method, as this may depend strongly on details of the implemen- 
tation. It is however reasonable to assume that all the methods considered in 
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Table 1 

Comparison of the implicit stability domains with full explicit stability for the 
schemes considered in this work and the traditional IMEX BDF2 and mCNAB 
methods. The results for the latter methods are taken from Reference [3]. 



Scheme 






a 


(16a) {k = 3) 






7r/2 


(16a) (k = 4) 






0.237r 


(16b) ip = 0, 


u = 1/3, k 


= 3) 


7r/4 


(16b) (/3 = 0, 


V = 1/3, k 


= 4) 


O.lSvr 


IMEX BDF2 






0.3l7r 


mCNAB (c = 


0) 







mCNAB (c = 


1/8) 




0.127r 


mCNAB (c = 


1/2) 




0.237r 



this study as well as the mCNAB and IMEX BDF2 schemes have comparable 
computational cost. The memory requirements may however be larger for the 
SSP-based schemes as data from more time levels need to be stored. 

We will also briefly discuss IMEX Rungc-Kutta methods. Ascher et al. [1] 
suggested that a two-stage, second-order IMEX Runge-Kutta scheme based 
on the trapezoidal rule could be an alternative to the CNAB method in 'most 
situations'. There are also schemes in this family with very good SSP pro- 
prieties. Analysis of IMEX RK2 schemes shows unfortunately that there are 
severe restrictions on the explicit stability domain, in particular for implicit 
eigenvalues with large magnitude [9]. Higher order IMEX RK methods have 
been proposed and analyzed by Ascher et al. [1], Pareschi and Russo [10], and 
Higueras [5], but these have not been considered in the present work. 
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